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STUDY TO DEVKWP GRADIOMETER TECHNIQUES 


I. INTRODUCTION 


There are three moving base gravity gradiometcrs currently under 
development. The instruments are being developed at Hughes Research 
Labs [Itef. l], the Bell Aerospace Corp. [2], and the Charles Stark 
Draper Lab. [:i, 1*), Tiie design goal for each of the sensors is 1 E6tv6*s, 

Furtheri.ioro, 0,1 Ebtvds (E) accuracy should be feasible from an orbiting 
gravity gradiometer [5]* The group of instruments includes sensors 
designed specifically to measure the gravity gradient, as well as sensors, 
which utilize existing accelerometers to provide a gradient estimate. 

The Hughes and Bell instruments rotate, thus modulating the information. 
This rotation transfers the gravity gradient signal to n higher fre- 
quency, quieter part of tlie spectrum, and can separate the signal from 
some sources of instrument bias. The Draper Lab sensor measures the 
gradient signal at zero frequency and uses a sophisticated flotation 
suspension system to isolate the sensing element from errors Induced 
by rotation and jitter. A system of at least three instruments of any 
one type is required in order to provide a complete gravity gradient 
tensor estimate. Often times, however, it is possible to extract all 
the information that is required from a single one of these instruments. 

The primary objectives of this paper, oriented toward the use of a 
workable gravity gradiometer as a sensing element In several applica- 
tions, are given below: 

1) To develop models for gravity gradient anomalies and gravity 
anomalies , 

2) To evaluate several methods of on-line instrument bias estima- 
tion, 

3) To determine the performance of a gradiometer in mapping the 
earth^s gravity field, 

4) To assess Inertial navigation systems augmented with a gravity 
gradiometer. 
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The work under this contract, and the work under a separate con-- 
tract (Goddard SPC NAS 5-2i960) to evaluate the performance 
of a geodesy mission with orbiting gradlometers have some overlap, 
especially in the groundwork areas of Instrument performance and capa- 
bility, and gravity models. As such, we hope it will be valuable to 
supply some results from work performed under that contract in this 
Final Report. 
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systems In which the ^jravity radiometer is an essential component, 
some model of the type of inputs is required. Several of the models 
are: 

1) tosserol harmonic models using Kaula's rule [7] , 

2) tesseral harmonic models using Allan's rule [8] , 

3) point mass and line mass models , 

4) experimentally determined second-order random process models. 


Once a give model is chosen, the gravitational potential, force, 
and gradient can be ascertained. The choice of which model is to be 
used in a given implementation depends on the dynamic range of interest 
for a given system (which is in turn dependent on the system's speed 
relative to the earth). This can vary substantially, from fixed base 
application, to ship speed, airplane speed, and finally, orbital speed. 
The various models also result in different degrees of complexity. 

This, too, must bo considered for any mechanization involving gravity 
grad ioinetcrs , 

This chapter contains the expressions for the gravity forces and 
gradients, and their correlation based on various gravity models. 


A. TESSERAL HARMONIC MODELS 


The most general expression of an arbitrary function over the 
surface of a sphere that satisfies the potential equation is in terms 
of tcsseral harmonics. The tcsscral harmonic expansion of the earth's 
gravitational potential is 


V(r,q),\) 
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whore (r,qp,^) is the position of a test point in terms of the spherical 

coordinates, radius, latitude, and longitude; m = GM ; = radius 

of the earth; P„ (x) are normalized associated Legendre polynomials, 

- ^ f m 

and C, and are normalized tesseral harmonic coefficients, which 

f m f m 

give the amount of each harmonic present; the perturbation potential 
due to a particular harmonic Is simply 


where 



( 2 ) 


and in writing (2), the phase Information is lost. 


Kstimates of the magnitudes of the given by Kaula [7] 

and Allan [Bl. Kaula *s familiar rule of thumb for is 


Im 
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(3) 


while Allan’s more complicated formula is 




Vl2,2 X 
^ ■ ■ ■ 


10 


-10' _„.(+3/2 


(0.93) 


(2( + 1) + 3 


(4) 


Although, different in form, both formulas agree quite closely over 
huge ranges of the subscript t, and with measured values of the 
tesseral harmonic coefficients^ 

With estimates of the magnitudes of the J. in (3) and (4), it 

{ m ' 

becomes possible to estimate the magnitude of the perturbation potential 

and hence of the force (i) and gradient (^) perturbations, since 
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If we further assume that the phases of the harmonics are random, 
it becomes possible to obtain closed-form approximate expressions for tlie 
variance of the perturbation forces, the variance of the perturbation 
gravity gradient, and the covariance of the perturbation forces and 
gradients due to the tesseral harmonics. 


The total variance of the force is computed for a representative 
force component (f^,) via Koula’s rule. It is given by 


2 2 2 2 2 
0^ ) y (f^ ) f (m sec ) 


( 6 ) 


where Is the radial force perturbation due to the ftli harmonic. 

The extra ( in (li) Is due to the fact there are f distinct harmonics 
(m =; 1 to t) for a particular f. 


fj. is obtained as the radial derivative 
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since 



where g is the gravitational acceleration at the surface of the earth. 

Squaring (7), inserting it into (6), using that pf (sin qp>) I = 1 

(in flvg 

by definition, and using Kaula^s rule for J. we obtain 

(m, 




.aj, . . X .0- E I iW 


(8a) 


so 


\ 10 


-5 


:(Wy ... 




(» 


2 V r 


(8b) 


An identical approach leads to closed-form expression for the grad- 
ients^ There 

rr \r/(m Cm 

2 'I 

where n - Using 
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eventually, using Kaula's rule. 
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Substitution of Allan's rule (4) into (D) gives a slightly more 
cumbersome series to sum. The results are 
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where 
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i.uxio"® r 3* 

<1 - x)^ 


= (o.,3^f . 


The results of Eqs. (8)^ (10), <ll), and (12) are plotted in Fig. 1. 
It shows the standard deviation of the force perturbations and gradient 
perturbations at various altitudes. Figure 1 also shows that as the test 
point approaches the earth*s surface, the standard deviation of the 
gradient using Kaula's rule blows up, a phenomena that physically does not 
occur. Although Kaula's rule Is in common use. it evidently does not 
attenuate the high frequency components rapidly enough to produce a finite 
standard deviation at the surface, Allan's rule does provide for a more 
rapid attenuation of the high frequency components, due to the presence 

p 

of the (0*93) term in (4), and does result in finite force and gradient 
variances , 

One useful feature of the closed-form solutions of the standard de- 
viations Is the ability to determine what amount of variance is due to 
different portions of the spectrum. A particularly simple example is 
obtained from the force components in (11). The analytical expression 
for due to harmonics above some degree k is easily found to be 

1.28 X 10 - X where again x = 0.93(R /r)^. Near the sur- 

face of the earth, then, the percent of variance contribution due to 
harmonic degree greater than or equal to k is given as 

k+2 

X 

1 - X 

2 
X 

1 - X 
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2 

where x ^ (0.9H). In order to account for only 50% of the noise, a 
model of fifth de^^ree 30 coefficients!) Is required, for 75% of 
tJie noise 10th (-'110 coefficients) degree, etc . The gradient components 
attenuate even less ropldly requiring an even higher order models This 
certainty proves that an improve<t tesseral harmonic model of the earth 
is not the way to go in order to Improve navigator accuracy due to the 
high di'grce of the model rct|uired to give even modest improvement. 

In order to compare results of other models to this model, the correl- 
ation co^'f f icien t s for P and g , and for P - P and g were 

rq^ rr qx\^ 

computed on the basis of the tessera 1 harmonic model with random phases. 
These correlation coefficients are given below near the earth’s surface 


and 



B, P01t>rr MASS AND LINE MASS MODE LS 

The point mass and line mass models approximate locai gravity 
anomalies much better than the tesseral harmonic model which has primary 
use in describing global features, 

A simplified analysis of n point disturbance begins with the assump* 
tion of a point mass with mass ~ Gm^, and a local 2-U x, y 

coordinate frame as is shown In Fig« 2« 


FIG. 2 POINT MASS MODEL M — ^ W 


t 


y 


- 9 - 


The (llsturbanco poteuLial duo to at a point (x^ y) is simply 

V - [l/v 0 js/'^ + . Expressions for the forces and gradients are 

easily obtained by taking the respective partial derivatives. 
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where r y/yp "•+ ^ . 

The CO r re la L ion b<' t wetMi the tour quant i t Ics In (13) can be found 
by assuming each of tliem is an independent measurement, and that the 
measuromcMits are obtained along a proscribed path In the x, y plane. 
For simplicity, the path chosen here is the path y constant, and x 
proceeds from -» to -i « at a constant rate (the path an airplane might 
follow). The information matrix is then given by 
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and the covariance matrix P Is simply 1 Normali/-ed off-diagonal 
elements of P give the correlation coefficients. A typical integra- 
tion is shown below 
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l-'roin ( 17 ) i\ is l’oun*i that V Is t'orrolatcd with f (aiui 

xy ^ ‘ _ X 

with nothiui; olso) with corro lat ItMi i'i>of flc it?tif 2 0*893^ and also 

that 1' - r is htiihlv corrolati’d with t (and with uoIIHuk olsr) 

XX yy ^ y 

with eor rol at. Hnt roiM'i lriont -4 -* O.H7f>, 

This lUMirly ptM’toot oorrolalton has particular si jrni f icanco for 
^ravily ^;radii>mot or system implement at ions , In most cases, the s^ravit y 
l»radiomeler is «se«l as a sensor to provide force perlurhai i ons due to 
Ki'uvily anomalit’s that the accelerometers cannot provide and that contain 
too miicli htith freqiieticy information for an a pr+orl earth K*’avity model 
to appVi>x imat Insleatl of proceedlntc alon^ the path of tntej^rnt tn^t the 
i;radietit Information to Rive the lorees (which leads to problems if t ht' 
grad ii>me ter measurement contains a Mas), the correliition of the forces 




with the gradients gives a way of proceeding directly from the gradients 
to the forces without augmenting the states of the system. In fact, 
from (16) we obtain 



<=^xy> 
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where y is a correlatioti height or a mean height of the vehicle above 
the disturbances. A simple gain adiustment is all that is needed then 
to provide the force perturbations from the gradient information. 


D. LINE MASS MODELS 

The line mass model is another possible alternate to locally model- 
ing gravitational perturbations. To determine if line mass models or 
point mass models result In better performance will require surveying 
data to discover if local gravity perturbations are more accurately 
modeled with lino masses or point masses. Surprlsi ngly ^ the use of 
point masses vs line masses produces relatively minor differences. The 
line mass model is shown in Fig. 3* 
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(19) 



The correlation between the various elements is obtained by the 
procedure described In (14) and (15), The identical substitutions are 
made here as in the point mass models the InteRrattons carried out^ and 


the rt'sults are that is correlated with (and with nothing else) 

with correlation coefficient \/j2 - 0,71 and that T - P is 

^ XX yy 


correlated with 
- \/j2 = -0.71. 


f 

y 


(and with nothing else) with correlation coefficient 


Although the correlation is not as strong as with the 


point mass models^ it is still quite high. Perhaps it is appropriate 
in the final model , to actually use a correlation coefficient between 
the values given by line mass models and point mass models. Again, 
scale heights can be derived to relate directly the gradient perturba- 
tions with the force perturbations. 
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An obvious advantage of these correlation models over dynamic models 
of the earth's gravity gradient perturbations ^ deflections of the vertical, 
and perturbation force components^ is the simplicity, A dynamic model 
involves earth surveying to ’’fit** parameters In the model, which may 
themselves vary substantially over different regions of the earth. 

An added advantage is that the correlation models do not involve *'aug" 
menting the state'* of the system, as do the dynamic models. 


Chapter III 
BIAS ESTIMATION 


Regardless of the eventual use of the gravity gradiometer, some type 
of signal processing will be required to minimize the effects of one of 
the largest gravity gradiometer errors: an unknown instrument bias. 

As will be discussed in Chapter V, a gradiometer bias has the same 
effect as a gyroscope drift when used in an Inertial navigator, l.e., 
they both produce unbounded position errors with Increasing time. When 
the gradiometer is used as the sensor to perform geodesy experiments 
from orbit, the gradiometer bias produces a bias in the estimate of the 
dominant terra of the earth's gravity field. 

This Chapter contains a study of several attempts to eliminate the 
gravity gradiometer bias error. 

One greatly simplified model of the bias estimation problem is 
given below. A single (rotating) gradiometer supplies two pieces of 
Information of the gravity gradient at a point, the difference of two 
principal elements of the gravity gradient tensor and one diagonal com- 
ponent . For example , 


*11 - *22 


= X 


12 


( 20 ) 


A set of three such gradiometers pointing in orthogonal directions 
along with Laplace's constraint equation (that the sum of the diagonal 
elements equal zero) is sufficient to determine the gravity gradient at 
a point; i.e., the problem below: 
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= 0 


( 21 ) 


is invertible^ and hence solvable for x in terms of z« (It is not 

T —1 T 

invertible without the constraint equation fact, S = (H H) Hz, 

Even if random noise v, with covariance r, is inserted into each of 

i i 

the first six equations in (21), it is possible to obtain a least - 

T —1 —1 T —1 

squares estimate for x in terms of z. In fact, x = (H R H) HR z 

T —1 —1 

and the variance of the estimate error P = (H R H) . Since Is 

an "exact" measurement thoimh, eltner Infinitesimal values for the 
corresponding element of R must be used and a limit process applies 
or else the matrix inversion lemma [10] may be used. The matrix inver- 
sion lemma can yield the Improved covariance matrix P from the covar- 
iance matrix P'. P* Is the easily computed covariance matrix which 
does not Include the beneficial effects of the exact constraint equation. 
The formula relating P to P’ Is 

P = P*H^(HP’i?’r^HP* 

where H reflects the constraint equation structure that 

"'22 *3 = 0 = ^ 

SO 

ii = [i 1 1 0 0 0], 
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The advantage of the matrix Inversion lemma procedure is that it 
completely avoids the limiting processes mentioned earlier. 

Now consider the problem of augmenting the state vector with con- 
•tant but unknown biases. 1 bias per gradiometer. Now 
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X + X + X =0 

11 22 33 


and although the sum of the biases b 4^ b b is obtainable (from 

1 iS 

z + z 4 z - z it is Impossible to obtain any other combination of 
the biases (in particular, the individual biases) and, therefore, it is 
impossible to solve for x. 


The only useful results from (23) are obtained if there is some 
a priori estimate of the biases. In this case, the usual Kalman filter 
techniques will allow for estimates of x to be obtained from noisy 
measurements, however the conditioning of the problem is of course de- 
pendent on the initial accuracy of b. If I is the initial informa- 

-1 ^ 

tion cm the augmented state (P^ ” ^0 then the update equations are 
simply 


+ (I^) H* R z 
T “1 

1^ = 4 R 


( 23 ) 


where is the measurement matrix for the augmented state vector. 
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and R the covariance of the measurement accuracy. 

Tiu' advantage of the initial Information, 1 , is evident in the 

T -1 ^ 

second equation of (23). Whereas H* K H* is not invertible, the sum 
T -1 

4 H* R H’ is invertible and hence the Inverses appearing in the estini- 
ate update equation are valid. 


More sophisticated schemes for estimating the biases make use of 
dynamic models of the system^ as opposed to the model in which both the 
gradient components and the biases were constants. Of course, since an 
a priori model of either the measurement or the state presumes additional 
Information Is available, it should be expected that better results could 
be obtained at the expense of additional complexity. 


Two additional bias estimation schemes will be discussed^ The first 
makes use of the fact that the bias is fixed in the instrument frame, while 
what the gradiometer measures will generally be fixed in an inertial frame. 
Rotating the spin axis of the gradiometer should then distinguish between 
the gradient components and the bias terms. The second method again 
distinguishes between the gradient components and the bias terms but now, 
the mechanism Is via a model of the gradient perturbations. This is where 
much of the material from Ch* I can be used. 


Beginning first with the rotation scheme, an expression for the 

gradient tensor is 
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where 
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10 0 
0 cos(uit) sin(uit) 


0 -sin(u)t) cos(u’t) 


and u> Is the rotation frequency. 


Then a single gradlometer . capable of measuring (say, without loss of 

generality) the T, , — P-,, and r,„ components would output 
XX 22 X2 

2 2 
= Tjj - (Pgj CO® wt + 2 P 23 am wt cos wt + sin uit) 


^2 = ^12 ’"is wt + bg + Vg 


where b and b are biases, and v and v are random noise. Now 
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( 24 ) 
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Again we have the constraint equation, ^ ^^33 Note 

that for t =. 0 , the t i me- varying H la (24) agrees with the H in 
(23), 

The standard test for observability can be performed on the matrix H 
in (24) to determine the maximal rank of the matrix 


0 



It turns out that the maximal rank of C is 7, so one mode is still 
not observable* Additional algebra yields the result that the mode 

is not observable, but that 2/3 b^^ + is observable. What 

this says is that it is still not possible to distinguish between and 

b* This result is easily explainable since the gradiometer spin axis 
never has a vertical component with the single axis rotation scheme w'p 
have considered here. However, if a uwo axis rotation scheme where the 
spin axis of the gradiometer spans all directions, all the parameters are 
observable, and hence the gradient components and Individual biases are 
obtainable from z. 

A problem with the continuous rotation schemes yet to be discussed 
is the introduction of a gradient field due to the kinematics of the rota- 
tion itself, The vector formula for this induced gradient field is 

" = [w x(tu X r)J = UKO - lu’ . (25) 

ui 

2 

In two dimensions the gradient field reduces to simply to . If the spin 

axis rotates as slowly as 2?r/50 sec, the required accuracy of this 

rotation rate to obtain full potential of the fixed instrument accuracy 
-9 

is 1 part in 10 I Since this accuracy is not feasible, the contin- 
uous rotation schemes are ruled out. 
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Ctlll, some use can be made of the preceding analysis. The fact that 
all the states and biases are observable with the continuous two axes 
rotation schema means that with a sufficient number of discrete measure^ 
ment points, the same information will be available. This can be accom- 
plished in several ways. One scheme is to ^'calibrate** the gravity gradio- 
meter at a fixed position, in several orientations to obtain the bias. If 
the bias is truly constant, the unit could then be used at new locations 
with the predetermined value of the bias. A second scheme is to simply 
repeat measurements with a movini^ base gravity gradiometer with the instru- 
ment oriented in different directions. For example, in a surveying mission, 
either with an airplane or satellite, it is possible to retrace a ground- 
track while making measurements in (three) different orientations. This 
determines the bias and states onr-llne, in which case the effects of bias 
drift could be minimized* 

The second scheme for distinguishing the gradient components from the 
biases is via a dynamic model of the gradient components. A simplifica- 
tion of a model to appear In a later chapter appears below. Again, letting 
X denote the gradient components, a linear dynamic model of x would be 
X = Fx + r*w where w is a random process. Since b = 0, there is again 
a mechanism for distinguishing x and b. Adjoining b to x 
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E(vv^) = R 
E(ww^) = Q . 
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T 

As bofor«, even ihoUK^ H II Ml in (26) not invert Iblo for use in the 
osllnuito updutr equation, the sum of the terms on the ri|cht-hand side of 
(26) is, and so eventually, It will be possible to distinguish x front 
b. More on the use of a model of the earth gravity field with applica- 
tloie- to initial iiuvigation will appeor later* 
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Clmptur IV 


GEODESY 


Si'verul expurimetiis have oraeri^eil which have the ciipubllit. y af very 
act'uralely nioa.surinK the ai'uvl tat lotml field of the earth. These n^^^t^lesy 
exper imt'iits which measure the hlKher harmonics of the earth Include (1) 
hlich-low mid low-low sutel I i t o-to-*siitel 1 1 to (S-S) trucking; (2) countor- 
orbi 1 1 ng (C ,0) d rag— f ree su^ o 1 1 i t es ; (3) u 1 1 Imeter measuroments f rom 

orbit; (-1) orbiting gravity gradlometer ntcasureineiits ; (0) others. 

The purpose of this chapter is to determine the accuriicy to which the 
higher harmonics of tlu' iMirth*s gravity field coulit be deduced with an 
orbiting gravity grad lometer , and to compare these results with some of 
t he pre V i ous 1 y ment loiied experiments , 

Till* starting point Is to expand the earth^s gravity fielil in a series 
of tosseral harmonies, ttecall from (1) the usual representation of 
the perturbation potentint us 


V(r.T,\) 


00 ( 11 / K ^ 

E E / /) V-*» ^>cs, 

f 2 m-0 


m 


cos m\ 


(27) 


-I S, sin m\l . 
tm 


It is possible to deduce niinlytic expressions for all the gradient compi>- 
iieiits, but for simplicity and since it will later be assumed that only 
r is measured, an expression for f will suffice here. (Actualtv 
there are six pieces of information available.) 


rr 


« ( 

57 J] ((4l)(f-t2) 
( 2 m 0 




(sin O’) [C cos m\ 
m t m 


(28) 


+ sin mV") . 

(m ^ 
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An Important point to be noted in (28) Is that although the entire 

gravity gradient tensor at a point cannot be reconstructed with Just 

r . It is possible to deduce values for all the tesseral harmonic co- 

rr 

efficients with Just (or in fact, with any single tensor component). 

In fact, the measurement i cp» X) + v is linear in the harmonic 

coefficients (with white noise v added) and so the whole theory of 
linear least-squares estimation theory can be applied. Letting the state 
vector X contain all tesseral harmonic coefficients, 


where 

T 

H 

((+1)(?+2)|— ^ j P^^(sin (p)cos ml. 

(f+l)(^+2)|^^— j P^^(sln <p) sin mX 
« 



* 

* 


(29) 
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®32* **• 
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®43’ ' 


®33' 

^^43 

T 

and surccsfiive rows of H 

are 

formed 

with increasing values of ^ and m* 


For satellite coverage of sufficiently large portions of the earth, the 
information matrix (inverse of the covariance matrix) can be approximated 
as 

. V' Hj^Cr,(p,Xjj) H <r,(p,X^) 

1 = 2 ^ a ( 30 ) 

K R 


since R is a scalar; (30) makes use of the fact that the satellite is in 
a polar orbit, and hence, latitude coverage is continuous, while longitude 
coverage is discrete . 


~ 2 &~ 




(ui)(t+2)(r+i)(r+2)|-^ j ) 


\2/ It 


P^n,(sin q,) P,,„,(8ln q,) 


f+( ' 


( ' nv' 


I cos m X, I I cos m* \ J 
[sin m X. j I sin 

RAcp 


(:n) 


Unless m - m* , then with sufficient coverage = 0, so harmonics with 

<lLfferent m decouple and the only remaining problem is 


2 H vt+tl ( 2 ) 


1 w 

u 


where T is the correlation time of the measurement noise in units of 
radians of orbit. Rewriting (32), we have 

2 

ii 2RT ^ r 

Unless I - t* the integral in (33) is zero. When the value 

of the integral is 


vt+t* ( co?»^nix| 

) jsln^raXj 


P. {sln<p)P,, (sliKp) cos cp <tp 
-till V m 


(33) 


« - »0m> 


SO I Is diagonal and the elements are 

X J 


^ a+i)^a+2)^ 

If Ir 

2RT ^ r r 


(*r i 


« - *0«> 


( 34 ) 


2 2 

where the average value of sin and cos is taken to be The 
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has diagonal eletnents 


information matrix for the C. and S. 

Im 4,m 


a+i)®(^+2)^ 


4 


n 


/R 

I ® I # orbits 1 2 

\ r / jr 2RT 2 

Om 


where 


2 

n 



( 35 ) 


This gives the parlance ' for a particular or of 




|m 


j 2«RT ' , 

fr ) 


f # orbits 


j 1 

* n^(t4l)(l+2) 


(36) 


Again, using Kaula's rule of thumb. 



14.14 


X 10“® . 


So, 


0-= ^ a 

tra. _ J 2it RT /r_\ 10 

J = »# orbit, * 14.14 


(37) 


Now 

n^ = 3000 E 

T - 10 sec = 0.0116 rad 

R = 0.1 E 

^ orbit = 3 months = 1200 orbits. 

Figure 4 shows a comparison of various three-month geodesy missions, 
including the current knowledge, counter orbiting satellites (C.O.), 
gravity gradlometers (G.G. ) at different altitudes, and satellite- to- 
satelllte (S-S) tracking. 
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Due to the fact that the gradioroeter yields a derivative measurenent , 
it tends to anpllfy the high frequency components of the earth's gravity 
field. This is evident in Fig. 4. The gradiometer outperforms the var- 
ious other experiments in the high harmonic range. 
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Chapter V 

KALMAN FILTERING FOR GRAVITY GRADIOMETER AUGMENTING 
INERTIAL NAVIGATION SYSTEMS 


A. INTRODUCTION 

Inertial navigation is based on the simple principle that position 
is given by the double integration of acceleration. Accelerometers do 
not measure accelerations themselves hut measure specific forces which 
are defined as all forces acting per unit mass with the exception of the 
gravity force. Hence, in order to obtain true accelerations, the specific 
force due to the gravity must be subtracted from the outputs of the accel- 
erometers, in other words, the gravity acceleration must be added: 

a = f + g (38) 

where 

0! s acceleration of the vehicle 

f = specific force 

g = gravity acceleration of the earth. 

Current inertial navigation systems use a "reference ellipsoid 
model” to compute the gravity acceleration of the earth, given the vehicle 
position. Although the reference ellipsoid model can well approximate 
the real gravity field of the earth, the difference is becoming a major 
error source of inertial navigation systems because of rapid hardware 
technology advances with accelerometers and gyroscopes [11], 

The difference between the actual gravity and the reference ellip- 
soid model may be expressed in terms of gravity anomaly (magnitude) 
and deflections of the vertical (angular deviation), [l2]. 
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From measurements taken at 12,5 nm intervals across the 35th parallel 
in the United States, the standard deviation and the correlation distance 
of the vertical deflection were determined as 5.2 arcseconds and 2S.lnml, 
respectively [Gelb, A. 1974]. The worldwide vertical deflection ensemble 
is considered to have 8 arcsecond rms and 20 n ml correlation distance [13]. 

Compensation for the errors caused by the gravity deflection and anom-* 
aly is one of the principal applications of gravity gradiometers currently 
under development. Gravity gradiometers measure gravity gradients which 
are related to gravity acceleration by the following relation: 

H « r • V (39) 


where 

r = gravity gradient tensor of the earth 
v = velocity of the vehicle. 

Given the gravity gradients with the velocity and the initial values 
of the gravity, we can integrate (39) to obtain the gravity. However, as 
is well known, bias errors of the gravity gradiometers may produce unbounded 
position errors as time increases. This fact is easily seen by the follow- 
ing simple example Illustrated in Fig. 5. Assiiming that only error source 
is the bias of the gr diometer, the linearized error propagation equation 
for a single horizontal channel with constant speed may be given by 


Ax a ilv 
Av =Ag^ 

Ag as vAF + r Av 

X XX XX 


(40a) 

(40b) 

(40c) 


where 

Ax. Av and Ag = estimation errors of position, velocity and 
gravity disturbance, respectively 
b = AF 5s bias error of the gradiometer. 

XX 
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For the spherical earth, the value of is 
- ^ ■ -1400 E. where R is the radius of the 
for F , we can integrate (40) to yield 


constant and equal to 
earth. Substituting - 



Ax 


Av 


= ^ t + -^ cos((o t + 0) 

2 U' s ^ 

U) s 

s 


vh 

“ + A sin((i-' t + 0) 


Aw cos(w t + 0) 
s s 


(4 la) 


(4 lb) 
(4lc) 


whore A and 0 are constants determined by 

vb 

+ A sin 0 a AVjj = initial velocity error 

ui 

s 

Aw^ cos 0 =1 ~ Initial gravity error. 

Equation (41a) indicates that the position error due to the bias of the 
gradiometer becomes unbounded with increasing time. 

In order to overcome this difficulty, Heller [13] proposed a method of 
using gravity gradiometer as an external aid combined with a gravity deflec- 
tion model. He obtained a number of numerical results for a single horizon- 
tal channel, assuming velocity reference errors, accelerometer errors, and 
gradiometer errors. 

In the following sections, we try to obtain an analytical solution for 
Heller's mechanization, considering the gradiometer error as the only error 
source. Then, we extend his mechanization to estimate the bias error of the 
gradiometer (the bias error in this case means the difference of the means 
of the outputs of the gradiometer from the gravity deflection model). 


B. SaLUTION FOR A SINGLE HORIZONTAL (mANKEL 

Th>! mechanization considered here is the same as "Gradiometer-as-an- 
external aid” (GAEA) in Heller [13], see Fig. 6. The gravity obtained 
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by the reference ellipsoid is used for navigation computation. External 
velocity information is provided in feedback form in order to damp the 
Schuler oscillation. Gradlometer measurements are combined with a gravity 
deflection model to estimate the navigation errors via use of a Kalman fil- 
ter. The vehicle speed is asstmied to be constant. 

Among various statistical gravity deflection models [Ref. 1S|, we 
choose the second order Markov model which is suitable for Kalman filter 
iropleiacntation because the governing equations may be written in the form 
of linear differential equations given by 


= -Pel + 

gi» = -pg|‘ + pWp 

where 

I = vertical deflection 
I* = augmented state 
^ S3 Markov parameter 

W. = zero mean white noise with power spectral density q. . 

6 5 

^ and are given by 

o 2.146V 

4g^llms 


(42) 


(43) 


wh'tre D and | are the correlation distance and rms of the vertical 
mis 

deflection respectively. 

In this section, we consider that the gradlometer error is described 
by a zero mean white noise. Then the measurement equation may be approxi- 
mated by 




= -Pg. + + vV 

I g 


(44a) 


36 - 


where 


r s output of gradiometer 

(r ) DC - ^ computed by using the reference ellipsoid model 

XX R£ XX 

B zero mean white noise of the gradiometer with power spectral 
density V^. 

V is given by 

D 

V„ = TAf^ (44b) 

B g 


where 

T = averaging time 

AFg rms of gradiometer error. 

For simplicity, and to make clear the effect of the gravity deflection, we 
assume that the only position error source is due to gravity deflection. 
Then, the error propagation equation for a single horizontal channel may 
be written as 


Ax = Av 

• 2 

Av a -W„Ax - 2ftiJ Av + g, 
8 ^8 


(45) 


where 

^ SB external velocity damping coefficient 
A R 

u n B Schuler angular frequency. 

S K 

As is clearly seen, the gravity deflection, |, and its augmented state |* 
are observable by the gradiometer measurement, but the navigation errors in 
position ^ix) and velocity CAv) are not observable. Hence, we can con- 

struct Kalman filter for | and |* and have the best estimates | and 

/\ 

g* given by 


= -Pg? + g|* -f K^(Z + vpgt - v3g?) 
g|*= -Pgi' + K^|(Z + vPg| - v3g|') 


(46) 
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Where K, and K., iirv Ktilnum filter izalns, 

Ni»w, the chilmution error In Rravity dlsturbnncv, RF.-gf'r drives 
the error equiit Ion (45), rather than the full K>^!ivity disturbance, gf,. 

The Kalman filter for the gravity deflection and Its augmented state 

may bo found as follows. First, the transfer function from the process 

noise w. to tlie measurement v. may be found as 
F. 


A 



S3 




(H + f^) 


2 


(47) 


Thon^ tho Hyinmeirlc root ohjirnclor 1st ic equation may ho written by 


or 


1 


(s -f v“V (-H + 0)^ 

K 


» 0 


1 + 


s 2 

;7 


(-s") 


(S + 1) 


(-S* + 1)^ 


= 0 


(4Ha) 


(48b) 


where 


F \ s 
s 

A (3 Pf. 12.140 

2v J V J nvi 

X g V 


rms 



* 


(4Kc) 

(48d) 


The symmetric root locus with paramoler a itiay bo readily drawn as shown 
in FIk* '7- this case, wo can solve the characteristic equations (48) 

and obtain the roots 


H 


H 


I 


I + a 


wo also find the steady-stale Kalman niter gains given by 


(49) 
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( 50 ) 


K 

K 


i 

I 


t 


s 0 



0 


« 


The cov.triunce of the estimation error of the gravity deflection is 

given by 





(51) 


As the iiceuracy of the gradiometer improves (AF 0 or a -> <») , the 

* ® 

error covariance P decreases monotonically (Pig. 8). For example, for 

f, r 

E =8 arcsec, D = 20 nrai, and v = 100 knots, we have from (43) 
rms 


P = 10.73 (hr"S 
q = 2.91 (n. rai^ hr ^) , 

K 


For T = 10 sec, AFg = IE , we have from (44b) 


V = 4.67 X 10~^(hr“^) . 

g 


Hence, from (48d) and (51), wo have 


a := 134 
P*^ = .0148, 

The rms of the estimation error in the vertical deflection decreases to 
about OHO tenth of that without the gradiometer. However, as a -»oo, 
one characteristic root approaches the origin and the other approaches 
infinity. This fact Indicates that the gravity deflection is obtained by 
integration of the measurement, ignoring the gravity deflection model. 
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MODAL DESTABIU2ATI0N 


0 €'& = 2.0 
A €/e = 0.1 
© €/p = 0.02 


without bias error 



FIG. 8 XCMRMALIZED COVARIANCE OF VERTICAL DEFLECTION ESTIMATICK^ ERRCR 






.N 

Substituting the ostimaition error of the gravity deflection |-^ for 
r in (dSb) and conducting tedious calculation, wo have expressions for 
the covariance of errors in position (P ) and velocity (P ) given by 

XX vv 


^ P 

A XX 


2(^+a^-l) + 4t;® yi+a* u)*^+4(; (1+a® )td* + yi+a® 'j 


XX (P > 
XX , 


P 

*^vv ■ <P ) 
vv 


2/4, 4,3 

* 


+ 4f;^u'** 4 4tcj* 4 1)A* 

fw S 

(52a) 

TFT ^ 

(52b> 


0 


a + 1) A*' 


where 


(p ) 

XX 0/ = covariances of errors in position and velocity, without 

. . I gradioDieter measurement 

' vv'o 


£s^ = u)*^ + 4(;,^+a® + (4f,®+4a^+2)co*® + 4tJ\+s!^ ui* + 1 


+ 2^i^* + 1) 

U)* = W /p . 

s s' • 


Among numerical examples shown in Figs, 9, the case with the vehicle veloc- 
ity 1000 knots is particularly interesting because when the root mean square 
values of the gravity gradiometer error are around 10 E, the covariances of 
position and velocity errors become larger than those without gradlometors. 
Ibis fact suggests studying the power spectral density of the estimation 
error of the gravity deflection which is given by 


(PSD) 




a 


*4 


U' 


2(l+2a*^)(/^ 


+ 1 


(53) 
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FIG. 9b NORMALIZED RMS VAIAJES OF ESTIMATI(»; ERRORS IN 
POSITIW, VELOCITY, AND VERTICAL DEFIECTION. 
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where 


A u) 

U' • 


Examining (53) and Fig* 10a we find that near zero frequency the power spec- 
tral density of the estimation error of the gravity deflection is always 
larger (up to a factor of 2) than q^. This means that even if the Kalman 
filter gives a smaller covariance as shown in Fig* 8^ it does not give im^ 
proved information near zero frequency, but worse. Hence, when the power 
spectral density of the gravity deflection error at Schuler frequency is 
larger than q^, the covariances of position and velocity errors are larg- 
€?r with gradiometers than without. 


C. BIAS ESTIMATE 

In this section, we make an attempt to extend Heller’s mechanization 
to estimate the bias error of the gradiometer, introducing the bias as an 
augmented state in the Kalman filter discussed in the previous section. 
The bias error in this case means the difference between the means of the 
ortputs of the gradlometers and the gravity deflection model* 

The measurement equation (44) may be rewritten as 

2 + Pgf,*+ vb + vVg (54) 

where b is the bias error of the gradiometer and the additional state 
equation is simply 


b = 0 (55) 

The full system consists of equations (42), (54) and (55)* Direct appli- 
cation of the Kalman filter theory to this problem falls, however* Since 
the bias Is an undlsturbable aid neutrally stable mode, the Kalman filter 
gain associated with this mode becomes zero after the initial uncortaintly 
disappears. This is a typical example of Kalman filter divergence. Many 
cures have been proposed for this difficulty such as restarting, minimum 
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FIG. 10a POWER SPECTRAL DENSITY OF ESTIMATION 


u (rad sec } 


tOR OF VERTICAL DEFIfCTION 









vitrinnce observers with eigenvalue constraints, added noise, pole-shifting 
and d<'Stabill7.atlon which are discussed In Bryson [14], Here, we use 
the modal destabilization method which Is based on the fact that the steady- 
stale Kalman filter for a system with an (undisturbed) unstable mode Is 
stable. As the amount of destabilization Increases, the absolute value of 
the eigenvalue associated with the undisturbed model increases from zero. 
However, the covariance of the estimation errors Increases too. Some num- 
erical results are shown In Fig. 8 and Table 1. The power spectral density 
of the estimation error of the gravity deflection was computed numerically 
and shown in Fig. lOb. Although there is still a hump higher than the orig- 
inal value q^, the power spectral density of 
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Table 1 


SOME NUMERICAL RESULTS 
MCmL 


OF SUBOPTIMAL FILTER OBTAINED BV 
DESTABILIZATION METHOD 


a 

e/p 

EIGENVALUES OF SUBOPTIMAL 
FILTER MEASURED IN UNITS OF 3 

0/4 ) 


0.001 

-200. 0 

-0.004, 

-0.003 

0.0243 

100 

0.01 

-200.0 

-0*125 ± 

jO.012 

0.073 


0.1 

-200.0 

-0.102 ± 

jO.084 

0.451 


2.0 

-20.0 

-0.527, 

-1.694 

0.979 


0.001 

-20.0 

-0.050, 

-0.004 

0.180 

10 

0.01 

-20.0 

-0.042, 

-0.028 

0.218 


0.1 

-20.0 

-0.124 ± 

JO. 098 

0.526 


2.0 

-21.9 

-0.356, 

-1.650 

0. 985 1 
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Chapter VI 

ESTIMATION OF DISTRIBUTED MASS DENSITY 


We have been trying a completely different way of processing gradio- 
meter measurement from that discussed in the previous chapter. Instead 
of equation (39) , we use the gradiometer measurement to estimate the mass 
density distribution of the earth, then compute the gravity from the den- 
sity distribution. This approach not only brings up a very interesting 
problem, i.e., filtering of a distributed system, but also seems to have the 
following practical advantages; 

(a) The bias of the gradiometer does not cause an unbounded position 
error but, at most, excites the Schuler oscillation. This is 
because both the gravity and the gravity gradients are obtained 
by spatial integration of the mass density and time-integration 
of hhe gradiometer measurement is not needed. 

(b) Since the statistical model of the mass density distribution is 
needed only at the boundary, the effect of the model error is 
small. Furthermore, there is little difficulty in extending this 
approach to the actual inhomogeneous earth density field. 

On the other hand, obviously this method requires large computer 
capacity. However, since the correlation distance of the gravity deflec- 
tion is about twenty nautical miles, we do not have to estimate the den- 
sity distribution over large areas for Inertial navigation purposes. This 
fact may relax the requirement of the computer capacity. 

So far, we have derived the partial differential equation for the 
density distribution and obtained the filtering algorithm, Including dis- 
tributed Kalman filter gain. Numerical calculation is in progress. 
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Chapter VII 

CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER STUDV 


Several gravity models, including 

1) tesseral harmonic models 

2) point mass and line mass models 

3) socoiul order random process models 

have been studied* As a result of the fact that Inertial navigators are 
particularly sensitive to what happens locally, it could be proved rigor- 
ously that the tesseral harmonic model probably could never be implemented 
due to the large number of parameters needed to obtain a suitably accurate 
local gravity description, 

Thu second order random process model was studied in detail. Wo ob- 
tained an analytical solution for Heller's mechanization (GAEA), using a 
more simplified problem formulation which still retains the part essential 
for gradlometor study, 1’he solution shows that the covariances of the 
errors in position and velocity are larger than those without gradiometer 
when Schuler frequency falls within the bamlv/ldth of the estimation error 
in the vertical deflection. 

We extended Heller's Kalman filter to estimate the bias error of the 
gradiometer, using the modal destabilization method to avoid Kalman filter 
divergence. 

Point mass and line mass models still should be considered, duo to 
the simplicity with which thev could be implemented. Future work will 
contain numerical results based on these models. 

Recently, another method of obtaining gravity perturbations has come 
under consideration, the estimation of distributed mass density. At this 
point, the necessary analytical work has been completed, and future work 
will include numerical results. 

The problcmt of bias estimation of the gravity gradiometer has now 
been studied in some detail. After several schemes were considered, the 
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most favorable methods to emerge are 

1) If possible, retrace a given groundtrack with the gradiometer 
in three different orientations. 

2) Calibrate the gravity gradiometer bias at a fixed location in 
many different orientations to obtain the bias accurately before 
the instrument Is used as a system component. 

Finally, the results of using a gravity gradiometer to perform a 
geodesy mission accurately are compared with competing schemes, A low 
orbiting gradiometer appears to be the most effective way of obtaining 
tlie tesseral harmonics of the earth for order 40 and above, l.e., the hig^ 
frequency components. In the range below the 40th harmonic, the geodesy 
mission employing counter-orbiting drag free satellites is superior. 
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